Optimal protocols and optimal transport in stochastic thermodynamics 
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Thermodynamics of small systems has become an important field of statistical physics. They 
are driven out of equilibrium by a control, and the question is naturally posed how such a con- 
trol can be optimized. We show that optimization problems in small system thermodynamics are 
solved by (deterministic) optimal transport, for which very efficient numerical methods have been 
developed, and of which there are applications in Cosmology, fluid mechanics, logistics, and many 
other fields. We show, in particular, that minimizing expected heat released or work done during 
a non-equilibrium transition in finite time is solved by Burgers equation of Cosmology and mass 
transport by the Burgers velocity field. Our contribution hence considerably extends the range of 
solvable optimization problems in small system thermodynamics. 
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The last two decades has seen a revolution in the un- 
derstanding of thermodynamics of small systenis driven 
out of equilibrium. Jarzynski's equality (JE) [1| relates 
an exponential average of the thermodynamic work W 
done on a system, driven from an initial equilibrium state 
to another final state, to the exponentiated free energy 
difference AF between these two states: 



Here and in the following /3 — l/ksT is the inverse tem- 
perature, ks the Boltzmann's constant and ^ • ;^ is an 
expectation over a non-equilibrium process, specified by 
a (time- and state-dependent) driving force or protocol. 
JE, and Crook's theorem 0, from which it follows, has 
been used to successfully determine binding free energies 
of single biomolecules through repeated pulling experi- 
ments [3], a feat which had previously been unimagin- 
able. For stochastic thermodynamics (the setting of this 
paper), such transient non-equilibrium fluctuation rela- 
tions are comprehensively reviewed in Q • A counter-part 
of the transient fluctuation relations are equally impor- 
tant steady-state fluctuation relations |5l-[llj. but these 
fall outside the scope of the present Letter where we con- 
sider only processes in a finite time interval. 

The transient non-equilibrium fiuctuation relations are 
identities; they hold irrespective of the protocol. Most 
quantities of interest however still depend on the pro- 
tocol, and can then be varied and optimized. A first 
step in this direction was taken by Schmiedl & Seifert 
who showed that when pulling a small system by opti- 
cal tweezers, (expected) heat released to the environment 
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and (expected) work done on the small system are mini- 
mized, not by naively smoothly pulling, but by protocols 
with discontinuities [12i | , a work which has generated con- 
siderable interest in the field For technical rea- 
sons, the analysis of Schmiedl & Seifert was limited to 
harmonic potentials. 

In this Letter we show how such optimization problems 
in stochastic thermodynamics (minimizing heat, work, 
the variance of the JE estimate of free energy differences) 
can be mapped to problems of (deterministic) optimal 
transport. The optimal control (for any of these cases) 
is determined by the solution of an auxiliary problem. 
When optimizing heat or work, this auxiliary problem is 
none other than the Burgers equation of fluid dynamics 
and cosmology, and mass transport by the Burgers field. 
Very efficient numerical methods have been developed to 
solve such problems, and these methods can be directly 
applied. Our contribution hence extends considerably 
the range of solvable optimization problems in stochastic 
thermodynamics. 

Stochastic thermodynamics and optimal protocols: We 
consider dynamics in the overdamped limit described by 
coupled Langevin equations: 

it = -\d^,V{$,,,t) + ^^w, , (2) 

with initial value ~ Xq, drift —d^^V and Wt a vector 
valued white noise with covariance {wtWf) — 5{t — t'). 
The mobility is and /3 the inverse temperature. For 
times t < to the potential is V{x,t) — Uo{x), and for 
times t > if is V{x,t) = U{{x). In the control interval 
[io,if] we allow the potential to be an explicit function 
of time V{x,t) — U{x,t) eventually discontinuous at the 
boundaries V{x,to) = poUo{x) + (1 — po)U{x,to) and 
V[x, tf) = pf U{x, tf) + (1 - pf ) Uf{x) with < Pi < 1, 
i = {o,f}. For single stochastic trajectories we define 
SW, the Jarzynski work |l|, and SQ, the heat released 



2 



into the heat bath, in as [l( 



5W= / dtV{^^,t) At , 



(3) 

(4) 



The difference SW — 5Q satisfies l*** law, i.e. is the inte- 
gral of an exact differential, if and only if the stochastic 
integral in 5Q is defined in the sense of Stratonovich. The 
Stratonovich integral is the limit of Riemann sums where 
products V{^,^,t)dwt are discretized according to to the 
mid-point prescription i.e. V{{(,^. + ^^.^^ )/2, ii)(wt;^j — 
uij.) for t G [ti, ti^i] and ti an arbitrary interpolation rule 
for t. Thus, the expression of the first law over [to, if] 



5W -5Q = V{£„^M)-V{^t^,to) 



(5) 



does not require pf = po = 1/2 for discontinuities in the 
time argument although the choice may appear otherwise 
appealing. 

The stochastic differential equations ^ lead to a 
(control-dependent) probability density m{x,t) evolving 
according the Fokker-Planck equation and the expecta- 
tion value of a local quantity Q is 



-< G{$t,t) ^= J m{x,t)g{x,t) 



(6) 



Straightforward application of Ito lemma (see e.g.[T3|) 
yields 



^dQ>^- j At < 4/^'^'^ ^ 



(7) 



for := ^ ■ dx + -^d^ the generator of the diffusion 
process with drift b. Given initial and final states, the 
minimal variance of the heat (or work) can be written 
as a KuUback-Leibler distance between a controlled and 
uncontrolled process, and this connection has been thor- 
oughly explored in the literature [H, [l^ . We will here 
be concerned with -< SQ -< 5W >- and exponentially 
weighted functionals of the heat or the work. 

Burgers equation in optimal stochastic control: We first 
focus on heat minimization. Following 20], we look for 
a function A{x,t) such that when evaluated along 







dt -< dtA 



u] 



(A + U) ^ 



(8) 



If such function can be found, the identity 

^ SQ^^^ A{^t^,t,)-A{^t^,ti) ^ , (9) 

holds true as -< {dt + £^^^*''^')A >- is the average of 
an exact stochastic differential. A sufficient condition for 
([5]) to be satisfied is the so-called dynamic programming 
equation (DPE) dtA + £L~^"^'(A + U) = which for 
any given value of U yields a linear, backwards in time 



evolution for A. The stationarity condition for DPE is 
obtained by taking the functional variation of ^ with 
respect to U. Introduce the (fictitious) potential R{x,t) 
corresponding to the state m{x,t) if it would have been 
in equilibrium i.e. R = logm. Then the variation of U 
yields the condition 



(10) 



which is satisfied independently of dxR if the potential is 

A-R , 



(11) 



where (j) is an arbitrary function of time alone. The op- 
timal control potential is therefore the solution of the 
coupled backwards- forwards equations 



ffl JiziA] A-\- R 



(12) 



dtm + dx ■ [^^^^m] = ^dlm . (13) 

ZT pT 

respectively obtained by plugging (fTTI) into the DPE 
and Fokker-Planck equations. We note that the Fokker- 
Planck equation has the property that if we split the drift 
into an equilibrium piece dxR and a remainder specified 
by the gradient of 



A + R 



(14) 



then it becomes the deterministic transport equation in 
the gradient of the remainder: 



dtm -\- -dx ■ [{dxip) m] = . 



(15) 



It is a perhaps surprising fact that using the definitions 
of R and ip and the Fokker-Planck equation (fT3|) reduces 
(ini) to simply 



dti^ + ^^=0 



(16) 



Equation (|T6)) is Burgers equation (for the velocity poten- 
tial), and equation ()15p is the equation of mass transport 
by the corresponding velocity field. These two equation 
are the first main result of this paper: we have reduced 
a complicated stochastic optimization problem to a clas- 
sical problem of optimal deterministic transport. In ad- 
dition, contrasting ([5]) with the expression of the work 
imposed by the first law, it is readily seen that work opti- 
mization brings about the same evolution equations (jl5l) , 
(1161) now complemented by the final boundary condition 

A{x,ti) = V{x,tt) = -[R{x,tf) +2i^{x,ti)] . (17) 

It is worthwhile remarking that the occurrence of final 
time constraints is a consequence of the backwards time 
evolution of the DPE and is a general feature of varia- 
tional principles in the presence of boundary cost terms 
[20j . We now discuss how this transport problem can be 
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solved, first if tlie initial and final states are given, and 
then if the initial state and the final control are given. 

Optimal heat between given initial and final states: The 
meaning of Burgers equation in (jl6p is somewhat peculiar 
in that it arises from a mixed forwards-backwards prob- 
lem. In other words, it is not reasonable to regularize 
(possible) shocks (in the future or in the past) by either 
adding -^-v&^ij: or —lyd^tp on the right hand side; equa- 
tion should make sense in both directions. On the 
other hand, without shocks the solutions of Burgers equa- 
tion are free-streaming motion, which we can specify by a 
inverse Lagrangean map Xq = X{ — {t{ — to)v{x{,t{) where 
the velocity (constant along streamlines) is v{xi,t{) = 
-dxfip{xf,tf). By mass conservation the inverse La- 
grangean map must satisfy the Monge- Ampere equation 



I det 



dx{ mo{xo) 



(18) 



where mo{x) = •m{x,to) is the initial state and mi{x) = 
m{x,tf) is the final state. In ID this equation is imme- 
diately solved in terms of the cumulative mass functions 
— ruf and = ruo- The inverse Lagrangean map 
is then determined by Mo{xo) = Aff(a;f). For higher 
dimensions we note that for free-streaming motion 



Xi 



tf~t 



*(a;f;tf,io)(J9) 



and (|18p becomes a partial differential equation in a 
scalar field vjf; 



det 



m[{x[) 



dxfdxf mo{da>,^) 



(20) 



Combining © with ([Til) the optimal released heat can 
be written as 

^SQ^^-^AS + 2< V'(lt„if)-V'(lt„,io) ^(21) 

where AS" = —(3 -< \nmi{^^J — ln77io(^j^) >- is the en- 
tropy change. Similarly, the minimal expected work is 

^ SW Vi^t^,tf) - V{^t^,t,) ^ + ^SQ^ , (22) 

provided R and ip satisfy ([T7)) . In both cases the dif- 
ference 2 ~< ip{^t^,tf) — ■!/'(^t^,io) >- which represents 
dissipated work, can also be written 



_^ Uu-n^.^^^^Qrr ^ ^^^^ 

t[ — to 



Equation (|23p means that the initial and final states 

x^^ 



can be specified by mass points {x'^\xl^\ . . . , Xg^''} and 



(1) ^(2) 
/ 



X) 



, x^^'*}, and a possible inverse Lagrangean 



{x 

map by a one-to-one assignment The inverse 

Lagrangean map solving (1^01) is then given by the assign- 
ment which minimizes the quadratic cost function (1231) . 
an approach which has been used with great success to 



reconstruct velocity fields in the early universe [21|, [22] ■ 
The interpretation of this quadratic cost function as dis- 
sipated work in stochastic thermodynamics is, up to our 
knowledge, new. 

Optimal heat with given final control: A setting which 
is closer to the problem of minimizing work discussed by 
Schmiedl and Seifert fl^ is when the final control V{x, tf) 
is specified, but not the final state. Let hence the initial 
potential be Uo{x) =|| x \\^ /2 + Uo and the initial state 
be mo{x) ^ exp{— /3 Uo{x)} and the final potential be 
Ui{x) = c \\ X - h \\^ /2 + U{ with c > and C7o , U{, 
arbitrary constants. By ([TT|) and ([T^ satisfies 



U4x,t) ^ -~[ijj{x,t) + R{x,t)] +(l){t) 



(24) 



The function can be, however, set to zero as the heat 
depends only upon the spatial gradient of J7*. Since for 
the heat there are no further conditions on tp and R, we 
can set U^{x,t{) = Uf{x) ~ V{x,tf). The problem can 
be then solved by a Gaussian Ansatz for the measure i.e. 



R{x,t) 



2 a? 



2/3 



In 



1 



(25) 



We can then use (|24p to write 4'{x,tf) in terms of ((25)) 
and Uf hence obtaining \I' by (fTO|) . Then, plugging ^ 
into ^ yields 



zTh 



& (7tf 



2T 



, (26) 



for T := tf — to and after straightforward algebra 

= ^/x,,. & a, = 1 + '-^{a,, - 1) , (27) 
for any t £ [to,t{]. Finally the optimal heat and drift are 



T 



T {at, - 1) 



T 



X t[x {atl - 1) + /XtJ 



Tat 



(28) 
(29) 



As expected the results do not depend on Uo ,U[. Whilst 
the state density ([25]) is continuous for all t G [to, tf], the 
optimal drift (I29p exhibits a discontinuity at t — to as 
discussed in [13| . 

Optimal work with given final control: Work optimiza- 
tion, as considered in exhibits more subtle features. 
The final condition (I17p together with (|24l) now yield 



ipix, tf) 



jl^ Pi) [Ui{x) + R{x,ti)] 



+ </'(if) .(30) 



Using the Gaussian Ansatz (j25p and proceeding as for 
the heat we find that dH]) still hold true but the 

final mean and variance are now given by 



cThpi 
PiTc + 2t 



ati 



K 



(31) 



with K = y/pfT{pfcT + 2T) + T'^ - r, and p; := (1 - 
pi) /{I — pi/2), i = {o, f}. As before, drift and density do 
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not depend upon (p nor Uo, U{. They, however, depend 
upon the shape of the discontinuities of the control V 
at the boundary. Note that for any c > 1, CTf,. is a 
decreasing function of T such that 1 > crtf > 1/ y/c. The 
corresponding expression of the optimal work is 



2rrf(l-gf,) 
/3T 

- AU , (32) 



4 (2 - po) I /3 CTt, 
T T Pi + 2 {I -at,) at, T 



where MJ = (1 - pf)f/f + pf 0(if ) - (1 - po) (/-(^o) - poUo 
can always be set to zero exploiting the arbitrariness of 
the function (f){t). It is straightforward to verify that the 
examples considered in [l2j are worked out for the case 
(po, pf ) = (1,0) and that as such they are a special case 
of the formulas given above. 

Optimizing the variance of the Jarzynski estimator: 
We now turn our attention to a different expectation 
value. The Jarzynski Equality (JE) is an equality in ex- 
pectationjp , but does not hold for a finite number of 
samples [2J]. Let there be N independent measurements 
of the work; then the free energy difference is estimated 
as AF = hi{j^ Ylf^i e"'^'^'), with a statistical er- 

ror determined by Var [e"^'^] /N. Moreover, expecta- 
tion and variance of a finite sampling will depend upon 
the details of the drift. It therefore makes sense to study 
the expectation value g\{x, t) =~< e~^^^ >-x,t, where we 
understand that the noise in the stochastic differential 
equation ([2]) is at inverse temperature /3, and that the 
initial state is in equilibrium at the same temperature. 
Using the approach of 9] g\ can be shown to satisfy for 
any given U a controlled diffusion equation which we can 
write for ^ ^"^S 3^ (note that go = 1 by defini- 
tion) as 

^[(2X^W]^^ = -£L(^-i)^-^lc/+ - II d^Ay^ f .(33) 

T 

The extremum condition for the drift gives 



(1 - 2A)Aa - i? 



(34) 



If we again, as in (jl4[) . split the drift into an equilibrium 
piece and a remainder 



9xV'A = -(l-2A)a, 



Ax + R 
2(1-A) 



(35) 



we obtain the generalized optimal transport equations 



dtm -f -dx ■ {mdxipx) = , 



5tV 



2t(1-2A) /3t(1-A) to/3t(A-1) 



(36) 
(37) 



These equations are not immediately solved, and deserve 
further study. 

In summary, we have shown how stochastic optimiza- 
tion problems are solved by the methods of optimal con- 
trol. The solution is built on an auxiliary problem of 
optimal transport. When minimizing heat or work of a 
small system this optimal transport is a classic of fluid 
mechanics and cosmology, namely Burgers equation. Be- 
tween any prescribed initial and final states, these prob- 
lems can be solved numerically with the Monge-Ampere- 
Kantorovich method, introduced to reconstruct velocity 
fields in the early Universe. Boundary cost contributions 
to the work, penalizing discontinuous controls and hence 
overcoming ambiguities in the definition of the free en- 
ergy, can be easily handled in the formalism in the form of 
Lagrange multipliers, and solved by fast-Legendre trans- 
forms methods for, at least, any convex potential. The 
direct connection between optimal transport and optimal 
protocols in small system thermodynamics was wholly 
unexpected, and is promising, as it applies to a whole 
wide class of related optimization problems. 
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